Here, we use Mathematica to investigate and graph the phase portrait of  the second example discussed in class for  the nonlinear systems.

({{x'= y}, {y'= x(1-(x^2))+y}})

Solve[{y==0, x(1-x^2)+y==0},{x,y}]

The equilibrium solutions are  (0,0), (-1,0) and (1,0)

J[x_,y_]=D[{y,x(1-x^2)+y},{{x,y}}]
Eigensystem[J[0,0]]
Eigensystem[J[-1,0]]
Eigensystem[J[-1,0]]

Based of the eigenvalues of J(0,0), the equilibrium solution (0,0) is a saddle. Based on the eigenvectors, the unstable manifold is tangent to  ((1/2) (-1+(Sqrt[5])),1) at (0,0) and the stable manifold is tangent to ((1/2) (-1-(Sqrt[5])),1) at (0,0). While based of the eigenvalues of J(±1,0), both equilibrium solutions (±1,0) are spiral sources.

graph1=VectorPlot[{y,x(1-x^2)+y},{x,-3,3},{y,-5,5},VectorPoints->20]

The vector plot seems to agree with our analytical results. Now we have to find solutions and graphs them to demonstrate the stable and unstable  manifold, and trajectories demonstrating the spiral trajectories. The difficulty is to pick the best initial data  to do this. But we can find and graph the solution both forward and backward in time.

For the stable manifold, I pick initial data on the vector ((1/2) (-1-(Sqrt[5])),1) and close to the origin. In positive time they should go to the origin and I choose just long enough time period to demonstrate their behavior. Note that below, the problem is solved for time -10 to 10, but I might not used the same time period  for each trajectory.

graph1=VectorPlot[{y,x(1-x^2)+y},{x,-3,3},{y,-5,5},VectorPoints->20] ;
gensol[x0_,y0_]:=NDSolve[{x'[t]==y[t],y'[t]==x[t](1-x[t]^2)+y[t],x[0]==x0,y[0]==y0},{x[t],y[t]},{t,-10,10}];
sol[1]=gensol[0.05(-1-Sqrt[5])/2,0.05];
sol[2]=gensol[ -0.05(-1-Sqrt[5])/2,-0.05];
graph2=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,2}]],{t,-10,2.25}, PlotRange->{{-3,3},{-5,5}},PlotStyle->{Thickness[0.01]}];
Show[{graph1,graph2},PlotRange->{{-3,3},{-5,5}},AxesLabel->{"x","y"}, Axes->True]

For the unstable manifold, I pick initial data on the vector ((1/2) (-1+(Sqrt[5])),1) and close to the origin. In positive time they should go to the origin and I choose just long enough time period to demonstrate their behavior. Note that below, the problem is solved for time -10 to 10, but I might not used the same time period  for each trajectory.

graph1=VectorPlot[{y,x(1-x^2)+y},{x,-3,3},{y,-5,5},VectorPoints->20] ;
gensol[x0_,y0_]:=NDSolve[{x'[t]==y[t],y'[t]==x[t](1-x[t]^2)+y[t],x[0]==x0,y[0]==y0},{x[t],y[t]},{t,-10,10}];
sol[1]=gensol[0.05(-1-Sqrt[5])/2,0.05];
sol[2]=gensol[ -0.05(-1-Sqrt[5])/2,-0.05];
sol[3]=gensol[0.1(-1+Sqrt[5])/2,0.1];
sol[4]=gensol[- 0.1(-1+Sqrt[5])/2, -0.1];
graph2=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,4}]],{t,-10,2.25}, PlotRange->{{-3,3},{-5,5}},PlotStyle->{Thickness[0.01]}];
Show[{graph1,graph2,graph2},PlotRange->{{-3,3},{-5,5}},AxesLabel->{"x","y"}, Axes->True]

But it looks like, I need longer positive time for the unstable manifolds and shorter negative time.  To do this I need to graph them separately.

graph1=VectorPlot[{y,x(1-x^2)+y},{x,-3,3},{y,-5,5},VectorPoints->20] ;
gensol[x0_,y0_]:=NDSolve[{x'[t]==y[t],y'[t]==x[t](1-x[t]^2)+y[t],x[0]==x0,y[0]==y0},{x[t],y[t]},{t,-10,10}];
sol[1]=gensol[0.05(-1-Sqrt[5])/2,0.05];
sol[2]=gensol[ -0.05(-1-Sqrt[5])/2,-0.05];
graph2=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,1,2}]],{t,-10,2.25}, PlotRange->{{-3,3},{-5,5}},PlotStyle->{Thickness[0.01]}];
sol[3]=gensol[0.1(-1+Sqrt[5])/2,0.1];
sol[4]=gensol[- 0.1(-1+Sqrt[5])/2, -0.1];
graph3=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,3,4}]],{t,-0.5,4}, PlotRange->{{-3,3},{-5,5}},PlotStyle->{Thickness[0.01]}];
Show[{graph1,graph2,graph3},PlotRange->{{-3,3},{-5,5}},AxesLabel->{"x","y"}, Axes->True]

Now, let’s add more solutions to show the spiral sources. To do this start close to (±1,0) and consider slight negative time and some positive time.

graph1=VectorPlot[{y,x(1-x^2)+y},{x,-3,3},{y,-5,5},VectorPoints->20] ;
gensol[x0_,y0_]:=NDSolve[{x'[t]==y[t],y'[t]==x[t](1-x[t]^2)+y[t],x[0]==x0,y[0]==y0},{x[t],y[t]},{t,-10,10}];
sol[1]=gensol[0.05(-1-Sqrt[5])/2,0.05];
sol[2]=gensol[ -0.05(-1-Sqrt[5])/2,-0.05];
graph2=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,1,2}]],{t,-10,2.25}, PlotRange->{{-3,3},{-5,5}},PlotStyle->{Thickness[0.01]}];
sol[3]=gensol[0.1(-1+Sqrt[5])/2,0.1];
sol[4]=gensol[- 0.1(-1+Sqrt[5])/2, -0.1];
graph3=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,3,4}]],{t,-0.5,4}, PlotRange->{{-3,3},{-5,5}},PlotStyle->{Thickness[0.01]}];
sol[5]=gensol[0.9, 0.1];
sol[6]=gensol[-0.9,- 0.1];
graph4=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,6}]],{t,-2,6}, PlotRange->{{-3,3},{-5,5}},PlotStyle->{Thickness[0.01]}];
Show[{graph1,graph2,graph3,graph4},PlotRange->{{-3,3},{-5,5}},AxesLabel->{"x","y"}, Axes->True]

A idea for changing the time interval for a solution is to read off a point on the solution and use it as the initial data. The point I chose on a spiral above was about (1.5, 2.9). I use it as the initial condition below.

graph1=VectorPlot[{y,x(1-x^2)+y},{x,-3,3},{y,-5,5},VectorPoints->20] ;
gensol[x0_,y0_]:=NDSolve[{x'[t]==y[t],y'[t]==x[t](1-x[t]^2)+y[t],x[0]==x0,y[0]==y0},{x[t],y[t]},{t,-10,10}];
sol[1]=gensol[1.5,2.9];
graph2=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,1}]],{t,-6,3}, PlotRange->{{-3,3},{-5,5}},PlotStyle->{Thickness[0.01]}];
Show[{graph1,graph2},PlotRange->{{-3,3},{-5,5}},AxesLabel->{"x","y"}, Axes->True]